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Abstract. We apply the transfer-matrix DMRG (TMRG) to a stochastic model, 
the Domany-Kinzel cellular automaton, which exhibits a non-equilibrium phase 
transition in the directed percolation universality class. Estimates for the 
stochastic time evolution, phase boundaries and critical exponents can be obtained 
with high precision. This is possible using only modest numerical effort since the 
thermodynamic limit can be taken analytically in our approach. We also point 
out further advantages of the TMRG over other numerical approaches, such as 
classical DMRG or Monte-Carlo simulations. 



1. Introduction 

The density-matrix-renormalisation-group (DMRG) is a powerful numerical tool which 
was developed by White ^ to investigate one-dimensional quantum chains. Since 
then it has been applied to several problems in quantum and classical physics |^ . 
DMRG studies of one-dimensional stochastic processes represent a relatively new field 
of interest. Such models can be described in terms of a non-symmetric stochastic 
Hamiltonian 10, The treatment of non-symmetric matrices within a DMRG 
algorithm [|[ a more challenging task than that of quantum systems. Carlon et 

al ^ have systematically studied the application of the "classical" DMRG to reaction- 
diffusion models. They showed that it is possible to obtain accurate results for non- 
equilibrium phase transitions and their critical properties. 

We propose here a different DMRG approach to study stochastic models, namely 
the transfer-matrix DMRG (TMRG). This method has first been applied to two- 
dimensional classical systems by Nishino Xiang et al [|o[ ^ used the TMRG 
to investigate the quantum transfer-matrix of one-dimensional quantum systems and 
introduced an accurate method for studying thermodynamic properties JT2[ . 
As an example we choose the Domany-Kinzel cellular automaton (DKCA) Q 
which exhibits a non-equilibrium phase transition in the universality class of directed 
percolation (DP) Its phase diagram and critical properties are well known 
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from previous studies, cf. p6[ for a review. We compare these literature data to 
our TMRG results to demonstrate their accuracy. Important advantages of our 
approach over DMRG and Monte-Carlo simulations (MCS) are pointed out, especially 
for systems like the DKCA for which the infinite time limit (t —>■ oo) and the 
thermodynamic limit {N — s- oo) do not commute. 

2. The Domany-Kinzel cellular automaton 

Cellular automata are algorithms that map one discrete configuration of a lattice to a 
new one in discrete time steps. Usually the map can be decomposed into local update 
rules. In the following we introduce the DKCA which is an example of a stochastic 
cellular automaton. 

2.1. The model 

The DKCA (p], |I| is defined on a periodically closed chain of N sites. Each site Si 
can take the values Si G {0, 1}. A configuration of the chain can be interpreted as a 
vector |s) = |si, . . . , sn) in a " Hilbert-space" . We speak of dead {si — 0) and active 
{si — 1) sites. The local update rules of a site Si{t + 1) are given by certain conditional 
transition probabilities 

p{s,{t + l)\s^^i{t)s,+i{t)) (1) 

which depend only on the neighbouring sites .Si±i{t). The model is controlled by three 
independent parameters 

po:=p(l|00), pi :=p(l|01) =p(l|10) and p2:=p(l|ll) (2) 

wherep(0|-) = 1— The DKCA is defined by po — 0. We are interested in the time 
evolution of an initial probability distribution P{t = 0) of lattice states. The phases of 
the model are characterised by the properties of the stationary distribution P{t = oo). 
The probability distribution P{t) can be regarded as a vector \P{t)) = J2s -^s{t)\s) 
where Ps{t) denotes the probability of a state |s) at time step t. 

2.2. Phases and Critical Phenomena 

The phase diagram differs, whether we consider a finite or infinite chain length N. 
For a finite system only one phase occurs for arbitrary pi,p2 < 1- Any initial state 
|P(0)) decays exponentially fast to the stationary state |0) :— |00 • ■ • O) which consists 
of dead sites only. |0) is referred to as an absorbing state of the model because the 
system cannot escape from |0) due to po — 0. However, the situation changes in the 
thermodynamic limit N —>■ oo. For sufficiently large pi andp2 an arbitrary initial state 
|P(0)) ^ |0) evolves into an unique stationary state |P(oo)) ^ |0). Thus we observe 
two phases, an absorbing and active one, separated by a sharp transition which falls 
into the DP universality class. 

Consider an arbitrary curve f{p) (0 < p < 1) in the phase diagram {pi,P2)- An order 
parameter for the transition from the active to the absorbing phase is given by the 
local density of active sites 

np{t) -.^ {l\n,\P{t)) (3) 
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with |l^ :— and the local number operator n^. In the absorbing phase np{t) 

vanishes in the limit t ^ oo, while in the active phase np{t) saturates at some 
stationary value np{oo). Close to the phase transition point pc a power law behaviour 

np{oo) {p - Pcf (4) 

is observed. In contrast to equilibrium systems without dynamical aspects, critical 
phenomena in non-equilibrium systems are usually characterised by two correlation 
lengths ill space and time respectively. Close to a phase transition point these 

correlation lengths diverge as 

^±-{p-Pcr-\ iw-iP-Pc)-"^^- (5) 

The DP universality class is already determined by the triple i^j^. ) ||l^. Two 
other critical exponents a and a are used in this work. Choosing p ^ pc, the dynamic 
evolution of np{t) scales as 

np(i)-t-", (6) 

and the response of np(oo) to an "external field" po as 

np,p,(oo)~(po)''/^ (7) 
a and a can be related to /?, ly^ and using scaling theory p!^ : 

f3 

a — — and a = + i'± — f3. (8) 
3. Transfer-matrix DMRG 

We first map the DKCA onto a two-dimensional classical system which can be solved 
using a transfer-matrix formulation. In the thermodynamic limit iV — > oo the 
properties of the model can be inferred from the knowledge of certain eigenvalues 
and -vectors of the transfer-matrix. The TMRG algorithm calculates this dominating 
part of the spectrum. 

3.1. Mapping 

The local update probabilities (|l|) of a site Si{t + 1) only depend on the neighbouring 
sites Si±i{t). Hence, we use a sub-lattice (parallel) update. An update time step 
t ^ t + 1 is divided into two sub-steps updating the sub-lattice of odd and even sites 
separately. Figure ^ depicts a realization of the update. The sites are connected by 
two state bonds. A bond is present (full line), if the predecessor Si±i{t) of a site 
Si(t -\- 1) is active, or not present (dotted line) otherwise. As shown in figure |l|, the 
bonds form a two-dimensional square lattice F. A configuration 7 G F of the lattice 
decomposes into vertices v which can be in eight different states. These are shown in 
figure H together with the probabilities assigned to them. 

Various properties of the stochastic system can now be expressed in terms of the 
statistical mechanics of a two-dimensional vertex model. For determining the critical 
behaviour we are mainly interested in the order parameter 

_ S,er"(7)p(7) 

Pil) — ^ denotes the probability of a certain lattice configuration 7 and 

corresponds to the statistical weight in a vertex model. 71(7) G {0,1} measures the 
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Figure 1. Sub-lattice parallel update of the DKCA. 



PO 






1 -Po 



Pi 



l-Pl 



Pi 



1-pi 



P2 



1 -P2 



Figure 2. Vertices of the classical model with corresponding probabilities. Full 
lines correspond to bonds. 



state of a certain lattice site, cf. figure Thus the calculation of np{t) requires taking 
the average (n) of a local operator n in a vertex model. Note that we have periodic 
boundary conditions in space direction N, but open ones in time direction t. 



3.2. Transfer-matrix formulation 

We briefly recall some basic facts of the transfer-matrix formulation of two-dimensional 
classical systems. The idea is to transform the sum over all configurations in (|^) into 
a product of transfer-matrices. The vertex model representing the dynamics of the 
DKCA has a checkerboard-like structure, cf. figure ^ Each dashed square tags a 
vertex and defines a local transfer matrix 
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The (vertical) transfer-matrix Tt of a certain time step t is defined by the matrix 
elements 

t 

{(72 ■ ■ ■ '72t\Tt\<j'2 ■ ■ •0'2t) =^ Y\. {'^2k-lCr2k\T\(72k-l<72k){<72k<72k+l\T\a2kO''2k+l)- (H) 

{Hk} k=l 
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In ( |11D the open boundary conditions in t-direction are considered by summing out 
fTi, CTi, o■!2l_^_l and a2t+i- Note that this definition is similar to the so-called quantum 
transfer matrix used in |Tl|, |l^ for quantum systems. In contrast to the horizontal 
transfer matrix in A^-direction, Tt is not a stochastic matrix. 

For the calculation of the order parameter another matrix Tt{n) is needed which 
contains one modified local transfer matrix (see figure |^) 



-(n) = 



/ 
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(12) 



Using Tf and Tt(ri) it can be easily shown that 

tr(T,(n)r,^/^-^) 



tr 



(13) 



Under the assumption that the highest eigenvalue Ag of Tf is not degenerate, the 
thermodynamic limit iV — > oo can be performed exactly 



np{t) 



tr(r.(n)T,^/^-i) {A^\T,H\A§) 



tr T, 



Ao 



(14) 



|A^) (|Aq )) denotes the right (left) eigenvector of the non-symmetric transfer matrix 
Tt corresponding to the eigenvalue Aq. 



3.3. TMRG algorithm 

The TMRG calculates Aq and the corresponding left and right eigenvector for 
an arbitrary matrix size t. This is done iteratively by numerical renormalisation 
techniques. We do not go into the details of the algorithm, which is basically the 
same used for the quantum transfer-matrix in ||ll[] . Two important modifications are 
necessary: 

• Open boundary conditions in time direction t are used in (^ij). These are 
necessary as the DKCA is not translationally invariant in time. In contrast, 
the corresponding Trotter dimension in quantum-TMRG is periodically closed 

[0. 

• The density matrix p^s — ^I'c \ -^o ){-^o \ used in quantum-TMRG is not symmetric 
where tro denotes the partial trace over the environmental block. Using p^s for 
the DKCA, we observed a poor convergence of the TMRG algorithm. Following 

we restricted ourselves to the symmetric density matrix 

p, = itre(|Ao^)(A«| + |A^)(Ao^|). (15) 

The convergence of the algorithm was improved remarkably. 

The eigenvectors |A^) and |Aq) of Tt were calculated by a power method. The 
diagonalisation of the density matrix was done by using Maple which can in principle 
handle arbitrary numerical precision. The data presented in the following section have 
been obtained by keeping to = 24 states. 
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4. Results 

4.I. Phase diagram 

In the absorbing phase np{t) decreases exponentially fast to zero, whereas rip (oo) ^ 
in the active phase. At the transition point p = Pc we expect a power law behaviour. 
Figure ^ (left) presents a double-logarithmic plot of rip along a curve f{p) in the 




Figure 3. The left figure shows the dynamic evolution of np{t) for different p 
on the SDP line. The curvature of the log/log-plot determines whether p belongs 
to the active or absorbing phase. The right figure depicts the ph a se diag ram 
calculated using TMRG (•) in comparison with data (o) from pa, n/A [isl DM. 



phase diagram (pi,p2)- As an example we used the site-directed percolation (SDP) 
line /sDp(p) = ip,p)- Pc is determined by adjusting p until the curvature of ripit) in 
the double-logarithmic plot vanishes. This was done for various lines fp2ip) = {p,P2) 
to calculate the complete transition hue, cf. figure^ (right). 



4-2. Order parameter 

In a second step we estimated np(oo) in the active phase by a 1/t finite-size scaling. 
Figure H depicts the results we obtained for the SDP line and the Hue fo{p) = (p, 0). 
Within the phases our data are in excellent agreement with MCS. Only close to the 
critical point we observe a small deviation. Here the 1 / t-scaling becomes less accurate 
due to critical slowing down. We expect that other methods, i.e. VBS or BST |20[ |2l|] , 
provide better results. 



4-3. Critical exponents 



The critical exponents a, (3 and cr determine the universality class, cf. section |2.2| . We 
exemplify the calculation of the exponents by investigating the SDP line. Figure H 
depicts for the order parameter (a) its time-dependence Upit), (b) the stationary value 
np(oo) and (c) the stationary value np^pg{oo) in the presence of a "field" po ^ withp 
close to Pc- The linear regression slope of a double- logarithmic plot (see insets) yields 
the particular exponent. Table |^ compares the results with other data obtained by 
series expansion and MCS. The accuracy of the exponents is similar to the DMRG 
calculations done by Carlon et al M. 
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Figure 4. TMRG results (•) of rip (oo) obtained from calculations across (a) the 
SDP-line and (b) the line /o(p) = (p, 0). MCS data (o) are plotted for comparison. 
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(this work) 





0.276486(8) 


0.27649(4) 


0.276(5) 




1.733847(6) 


1.73383(3) 






1.096854(4) 


1.09684(1) 




a 


0.159464(5) 


0.15946(2) 


0.156(5) 




0.10825(2) 


0.10825(2) 


0.108(5) 



Table 1. Estimates for the critical exponents obtained by series expansion and 
MCS in comparison with TMRGL Results for and can be obtained for our 
data using the scaling relations (H). 



5. Conclusions and Outlook 

In this letter we have proposed a new approach to stochastic models using TMRG. 
The accuracy of the results were shown by comparison to MCS or literature data. 
Critical exponents were estimated which deviate from the currently accepted values 
[ p2[ by less than 3%. In the following we focus on three important advantages of 
TMRG compared to traditional MCS or DMRG calculations. 

First, the exact treatment of the thermodynamic limit of the chain is the main aspect 
of the TMRG presented here. In contrast, MCS and DMRG can only be applied 
to finite chains. Thus a finite-size scaling has to be used carefully to estimate the 
thermodynamic limit N ^ oo. For transitions into absorbing states N oo does 
usually not commute with the infinite time limit i — > c». As shown in N —> oo 
can be performed exactly within the transfer-matrix formulation. 
Second, the TMRG approach facilitates a proper investigation of the dynamic 
evolution of the system (i.e. np{t)). This is a more difficult task using MCS. Many 
samples have to be taken into account to obtain a reasonable average. In contrast, the 
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Figure 5. Calculation of the critical exponents (a) a, (b) f3 and (c) ^. The 
insets are double-logarithmic plots of the particular data. The critical exponents 
are calculated by determining the linear regression slope of the log/log-plot. 



DMRG and TMRG approaches do not require random numbers and averaging like the 
MCS. The DMRG algorithm basically calculates the ground state of the stochastic 
Hamiltonian H which corresponds to the steady state (t — > oo) Dynamic properties 
are available by analysing the gap of H. Thus the second eigenvalue of H has to be 
computed. Critical properties are then extracted from the finite-size data of the gap 
[^. On the other hand, the TMRG offers direct access to the dynamics and phase 
transition of the model by calculating one eigenvalue only. 

The third advantage is a rather technical one. We obtained surprisingly accurate 
results by modest numerical efforts. In comparison, the classical DMRG algorithm 
needs more sophisticated numerical techniques to guarantee a sufficient convergence 
[^. We expect a further improvement of our results by involving various numerical 
refinements, i.e. a modified Arnoldi method 



21 



or a EST finite size scaling |2C 
The TMRG approach can be generalized to study other stochastic processes with 
arbitrary dynamics. Currently we investigate the pair contact process with diffusion 
(PCPD) and compare TMRG and DMRG results ||||. The universality class of the 
model is still unknown and an interesting field of research E^, EH] . 
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